Scoliosis Model

Autor

Caio Vallio

Data de Publicação

23 de dezembro de 2025

Resumo

Contexto: A escoliose idiopática do adolescente pode apresentar resposta clínica variável a intervenções conservadoras. Objetivo: Explorar associações entre características basais e melhora radiográfica em 6 meses. Métodos: Estudo observacional com análise exploratória e inferencial. O desfecho contínuo foi \(\Delta\) (maior curva em 6 meses − baseline, em graus). O desfecho binário (MCID) foi melhora definida como \(\Delta \le -5°\). Ajustamos um modelo linear para \(\Delta\) e um modelo logístico para o MCID, reportando estimativas com IC95% e verificações de adequação dos modelos.

0.1 Como reproduzir

  • O arquivo de dados é lido de data/modelagem final.xlsx (aba dados).
  • As análises dependem de pacotes R listados no chunk de setup; caso algum pacote esteja ausente, o documento interrompe a execução com uma mensagem explícita.

0.2 Objetivo e delineamento

Este é um relatório exploratório e inferencial: o objetivo é estimar associações ajustadas entre variáveis basais e a evolução radiográfica em 6 meses.

0.3 Definições de desfecho

  • Desfecho contínuo: \(\Delta\) = (maior curva em 6 meses − maior curva baseline), em graus. Valores mais negativos indicam maior melhora.
  • Desfecho binário (MCID): melhora clínica definida como \(\Delta \le -5°\) (redução de pelo menos 5 graus).

0.4 Plano de análise estatística

  • Descrição da amostra: estatísticas descritivas das variáveis basais e do desfecho.
  • Inferência (associações ajustadas):
    • Modelo linear para \(\Delta\) (efeitos como betas, IC95%).
    • Modelo logístico para MCID (efeitos como odds ratios, IC95%).
  • Adequação dos modelos: diagnósticos de assunções (colinearidade, resíduos, heteroscedasticidade quando aplicável) e medidas de qualidade de ajuste (por exemplo, \(R^2\)).
  • Forma funcional: para preditores contínuos, assume-se relação aproximadamente linear.

Por se tratar de análise exploratória, os resultados devem ser interpretados com cautela, com ênfase em magnitude/direção e incerteza (IC95%).

Código
## Setup (reprodutibilidade)
pkgs <- c(
    "tidyverse", "gtsummary",
    "readxl", "janitor",
    "performance", "qqplotr", "PupillometryR"
)
missing_pkgs <- pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_pkgs) > 0) {
    stop("Pacotes ausentes: ", paste(missing_pkgs, collapse = ", "))
}
invisible(lapply(pkgs, library, character.only = TRUE))

set.seed(123)
theme_set(theme_minimal())


# read data
df_raw <- readxl::read_excel("data/dataset_escoliose_00.xlsx", sheet = "dados", na = "")

# clean names
df <- df_raw |> janitor::clean_names()

# escoliometro - classificacao do local onde a escoliose esta mais severa
df <- df |>
    mutate(
        escoliometro_maior_10_graus = case_when(
            (escoliometro_cervical > 10) & (escoliometro_torarica > 10) & (escoliometro_lombar > 10) ~ "cervical_toracica_lombar",
            (escoliometro_cervical > 10) & (escoliometro_torarica > 10) ~ "cervical_toracica",
            (escoliometro_cervical > 10) & (escoliometro_lombar > 10) ~ "cervical_lombar",
            (escoliometro_torarica > 10) & (escoliometro_lombar > 10) ~ "toracica_lombar",
            (escoliometro_cervical > 10) ~ "cervical",
            (escoliometro_torarica > 10) ~ "toracica",
            (escoliometro_lombar > 10) ~ "lombar",
            TRUE ~ "nenhum"
        )
    )

# lenke - classificacao da curvatura da escoliose
df <- df |>
    mutate(
        lenke_curvatura = case_when(
            lenke %in% c(1, 2, 5) ~ "curva em c",
            lenke %in% c(3, 4, 6) ~ "curva em s",
            TRUE ~ NA_character_
        )
    )

# risser - classificacao do inicio e final do crescimento
df <- df |>
    mutate(
        risser_agrupado = case_when(
            risser %in% c(0, 1, 2) ~ "inicio do crescimento",
            risser %in% c(3, 4) ~ "final do crescimento",
            TRUE ~ NA_character_
        )
    )

# cifose - classificacao da cifose
df <- df |>
    mutate(
        cifose_categ = case_when(
            cifose_toracica < 10 ~ "hipocifose",
            cifose_toracica >= 10 & cifose_toracica < 40 ~ "normalidade",
            cifose_toracica >= 40 ~ "hipercifose",
            TRUE ~ NA_character_
        )
    )

# lordose - classificacao da lordose
df <- df |>
    mutate(
        lordose_categ = case_when(
            lordose_lombar < 45 ~ "hipolordose",
            lordose_lombar >= 45 & lordose_lombar < 50 ~ "normalidade",
            lordose_lombar >= 50 ~ "hiperlordose",
            TRUE ~ NA_character_
        )
    )

# definicao do desfecho - alvo
# cobb inicial - alvo
df <- df |>
    mutate(
        cobb_inicial_maior = pmax(
            cobb_toracico_proximal,
            cobb_toracica,
            cobb_lombar,
            na.rm = TRUE
        )
    ) |>
    mutate(
        regiao_cobb_inicial = case_when(
            cobb_inicial_maior == cobb_toracico_proximal ~ "toracico_proximal",
            cobb_inicial_maior == cobb_toracica ~ "toracica",
            cobb_inicial_maior == cobb_lombar ~ "lombar",
            TRUE ~ NA_character_
        )
    )

# categorizacao do desfecho - alvo
df <- df |>
    mutate(delta = maior_curva_6_meses - cobb_inicial_maior) |>
    mutate(delta_cat = if_else(delta <= -5, 1, 0)) |>
    mutate(
        delta_cat_f = factor(
            delta_cat,
            levels = c(0, 1),
            labels = c("Sem melhora (MCID)", "Melhora (MCID)")
        )
    )

# casting de variaveis (garantir tipos adequados no ajuste)
df <- df |>
    mutate(
        idade = as.double(idade),
        lenke = as.factor(lenke),
        risser = as.factor(risser),
        sexo = as.factor(sexo),
        flexibilidade = as.factor(flexibilidade),
        regiao_cobb_inicial = as.factor(regiao_cobb_inicial),
        escoliometro_maior_10_graus = as.factor(escoliometro_maior_10_graus),
        lenke_curvatura = as.factor(lenke_curvatura),
        risser_agrupado = as.factor(risser_agrupado),
        cifose_categ = as.factor(cifose_categ),
        lordose_categ = as.factor(lordose_categ)
    )

0.5 Dados faltantes e tamanho amostral

Código
missing_tbl <- df |>
    summarise(across(everything(), ~ sum(is.na(.)))) |>
    pivot_longer(everything(), names_to = "variavel", values_to = "n_missing") |>
    mutate(pct_missing = n_missing / nrow(df)) |>
    arrange(desc(pct_missing))

missing_tbl |>
    mutate(pct_missing = scales::percent(pct_missing, accuracy = 0.1)) |>
    print(n = 50)
# A tibble: 29 × 3
   variavel                    n_missing pct_missing
   <chr>                           <int> <chr>      
 1 escoliometro_cervical             262 97.8%      
 2 cobb_toracico_proximal            219 81.7%      
 3 cobb_lombar                        36 13.4%      
 4 escoliometro_lombar                27 10.1%      
 5 cobb_toracica                      23 8.6%       
 6 flexibilidade                      16 6.0%       
 7 escoliometro_torarica              13 4.9%       
 8 id                                  0 0.0%       
 9 idade                               0 0.0%       
10 sexo                                0 0.0%       
11 altura                              0 0.0%       
12 peso                                0 0.0%       
13 lenke                               0 0.0%       
14 risser                              0 0.0%       
15 cifose_toracica                     0 0.0%       
16 lordose_lombar                      0 0.0%       
17 dif_colete                          0 0.0%       
18 correcao_colete                     0 0.0%       
19 maior_curva_6_meses                 0 0.0%       
20 escoliometro_maior_10_graus         0 0.0%       
21 lenke_curvatura                     0 0.0%       
22 risser_agrupado                     0 0.0%       
23 cifose_categ                        0 0.0%       
24 lordose_categ                       0 0.0%       
25 cobb_inicial_maior                  0 0.0%       
26 regiao_cobb_inicial                 0 0.0%       
27 delta                               0 0.0%       
28 delta_cat                           0 0.0%       
29 delta_cat_f                         0 0.0%       
Código
# drop na in flexibilidade
df <- df |> drop_na(flexibilidade)

Observação: os modelos abaixo usam análise por caso completo (listwise deletion), que é o comportamento padrão do glm()/lm() quando há dados faltantes nas variáveis do modelo.

0.6 Codificação de variáveis categóricas (referências)

As variáveis categóricas entram como fatores. A categoria de referência (baseline) é o primeiro nível do fator (conforme levels()), a menos que seja explicitamente reordenado.

Código
categoricas <- c(
    "sexo",
    "lenke",
    "risser",
    "flexibilidade",
    "regiao_cobb_inicial",
    "escoliometro_maior_10_graus",
    "lenke_curvatura",
    "risser_agrupado",
    "cifose_categ",
    "lordose_categ"
)
map_dfr(categoricas, function(v) {
    tibble(
        variavel = v,
        niveis = paste(levels(df[[v]]), collapse = " | ")
    )
}) |>
    print(n = Inf)
# A tibble: 10 × 2
   variavel                    niveis                                      
   <chr>                       <chr>                                       
 1 sexo                        feminino | masculino                        
 2 lenke                       1 | 2 | 3 | 4 | 5 | 6                       
 3 risser                      0 | 1 | 2 | 3 | 4                           
 4 flexibilidade               flexivel | rigido                           
 5 regiao_cobb_inicial         lombar | toracica | toracico_proximal       
 6 escoliometro_maior_10_graus lombar | nenhum | toracica | toracica_lombar
 7 lenke_curvatura             curva em c | curva em s                     
 8 risser_agrupado             final do crescimento | inicio do crescimento
 9 cifose_categ                hipercifose | hipocifose | normalidade      
10 lordose_categ               hiperlordose | hipolordose | normalidade    




1 Análise descritiva

1.1 Variáveis de entrada do modelo

Variáveis coletadas na linha de base.

Código
df |>
    gtsummary::tbl_summary(
        include = c(
            # numericas
            idade,
            altura,
            peso,
            cifose_toracica,
            lordose_lombar,
            dif_colete,
            correcao_colete,

            # categoricas
            sexo,
            lenke,
            risser,
            flexibilidade,
            regiao_cobb_inicial,
            escoliometro_maior_10_graus,
            lenke_curvatura,
            risser_agrupado,
            cifose_categ,
            lordose_categ,
        ),
        type = list(
            idade ~ "continuous"
        ),
        statistic = list(
            gtsummary::all_continuous() ~ "{mean} ({sd})",
            gtsummary::all_categorical() ~ "{n} ({p}%)"
        )
    )
Characteristic N = 2521
idade 12.96 (1.58)
altura 1.58 (0.08)
peso 49 (27)
cifose_toracica 22 (11)
lordose_lombar 54 (12)
dif_colete -17.4 (6.4)
correcao_colete 48 (19)
sexo
    feminino 225 (89%)
    masculino 27 (11%)
lenke
    1 14 (5.6%)
    2 18 (7.1%)
    3 56 (22%)
    4 23 (9.1%)
    5 24 (9.5%)
    6 117 (46%)
risser
    0 54 (21%)
    1 43 (17%)
    2 65 (26%)
    3 49 (19%)
    4 41 (16%)
flexibilidade
    flexivel 158 (63%)
    rigido 94 (37%)
regiao_cobb_inicial
    lombar 137 (54%)
    toracica 112 (44%)
    toracico_proximal 3 (1.2%)
escoliometro_maior_10_graus
    lombar 88 (35%)
    nenhum 52 (21%)
    toracica 94 (37%)
    toracica_lombar 18 (7.1%)
lenke_curvatura
    curva em c 56 (22%)
    curva em s 196 (78%)
risser_agrupado
    final do crescimento 90 (36%)
    inicio do crescimento 162 (64%)
cifose_categ
    hipercifose 20 (7.9%)
    hipocifose 25 (9.9%)
    normalidade 207 (82%)
lordose_categ
    hiperlordose 169 (67%)
    hipolordose 50 (20%)
    normalidade 33 (13%)
1 Mean (SD); n (%)

1.2 Desfecho principal

  • Maior curva na linha de base e em 6 meses.
  • Delta: maior curva 6 meses - maior curva baseline.
  • Delta categórica: delta melhor em 5 graus ou mais (MCID).
Código
df |>
    gtsummary::tbl_summary(
        include = c(
            cobb_inicial_maior,
            maior_curva_6_meses,
            delta,
            delta_cat_f
        ),
        type = list(
            delta_cat_f ~ "categorical"
        ),
        statistic = list(
            gtsummary::all_continuous() ~ "{mean} ({sd})",
            gtsummary::all_categorical() ~ "{n} ({p}%)"
        )
    )
Characteristic N = 2521
cobb_inicial_maior 36.7 (5.6)
maior_curva_6_meses 32 (8)
delta -4.5 (5.5)
delta_cat_f
    Sem melhora (MCID) 135 (54%)
    Melhora (MCID) 117 (46%)
1 Mean (SD); n (%)

1.3 Distribuição do desfecho principal

Código
df |>
    ggplot(aes(x = "", y = delta)) +
    PupillometryR::geom_flat_violin(
        position = position_nudge(x = .2),
        # fill = "steelblue",
        alpha = 0.7
    ) +
    labs(
        title = "Distribution of the outcome (delta)",
        x = "",
        y = "delta (degrees)"
    ) +
    geom_point(position = position_jitter(w = .15)) +
    geom_boxplot(
        width = .25,
        alpha = 0.7,
        outlier.shape = NA
    ) +
    coord_flip() +
    theme(
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()
    )

Código
df |>
    ggplot(aes(sample = delta)) +
    qqplotr::stat_qq_band(alpha = 0.2) +
    qqplotr::stat_qq_line() +
    qqplotr::stat_qq_point() +
    labs(
        title = "Outcome QQ-plot (delta)",
        x = "Theoretical quantiles",
        y = "Sample quantiles"
    )




2 Modelagem

2.1 Modelo de regressão linear

Este modelo estima associações ajustadas com o desfecho contínuo \(\Delta\) (maior curva em 6 meses − baseline, em graus).

Código
model_lin <- lm(
    delta ~
        # numericas
        idade +
        altura +
        peso +
        cifose_toracica +
        lordose_lombar +
        dif_colete +
        correcao_colete +

        # categoricas
        sexo +
        lenke +
        risser +
        flexibilidade +
        regiao_cobb_inicial +
        escoliometro_maior_10_graus +
        lenke_curvatura +
        risser_agrupado +
        cifose_categ +
        lordose_categ,
    data = df
)

2.1.1 Parâmetros do modelo

Código
model_lin |>
    gtsummary::tbl_regression(conf.level = 0.95) |>
    # gtsummary::add_global_p() |>
    gtsummary::bold_p() |>
    gtsummary::bold_labels() |>
    gtsummary::italicize_levels()
Characteristic Beta 95% CI p-value
idade -0.11 -0.63, 0.41 0.7
altura 0.09 -9.3, 9.5 >0.9
peso -0.01 -0.04, 0.01 0.3
cifose_toracica -0.07 -0.17, 0.02 0.11
lordose_lombar 0.08 -0.01, 0.17 0.10
dif_colete 0.27 0.04, 0.50 0.023
correcao_colete -0.01 -0.08, 0.06 0.8
sexo


    feminino
    masculino -0.44 -2.6, 1.8 0.7
lenke


    1
    2 2.4 -1.0, 5.7 0.2
    3 5.0 2.2, 7.9 <0.001
    4 6.2 2.9, 9.4 <0.001
    5 1.3 -4.6, 7.2 0.7
    6 4.0 -1.6, 9.5 0.2
risser


    0
    1 -2.5 -4.5, -0.47 0.016
    2 -1.1 -3.0, 0.73 0.2
    3 -0.24 -2.2, 1.7 0.8
    4 0.13 -2.3, 2.6 >0.9
flexibilidade


    flexivel
    rigido 1.7 0.37, 3.0 0.012
regiao_cobb_inicial


    lombar
    toracica 0.57 -4.4, 5.5 0.8
    toracico_proximal -3.6 -11, 4.1 0.4
escoliometro_maior_10_graus


    lombar
    nenhum 0.39 -1.4, 2.2 0.7
    toracica 0.53 -1.4, 2.5 0.6
    toracica_lombar 2.6 -0.05, 5.2 0.054
lenke_curvatura


    curva em c
    curva em s


risser_agrupado


    final do crescimento
    inicio do crescimento


cifose_categ


    hipercifose
    hipocifose -3.0 -7.7, 1.8 0.2
    normalidade -2.6 -5.8, 0.72 0.13
lordose_categ


    hiperlordose
    hipolordose 2.4 -0.12, 4.9 0.062
    normalidade 1.2 -1.1, 3.4 0.3
Abbreviation: CI = Confidence Interval

2.1.2 Verificações de adequação do modelo

2.1.2.1 Colinearidade

Verificação de colinearidade entre as variáveis independentes (Variance Inflation Factor). Valores maiores que 10 indicam colinearidade entre as variáveis.

Código
performance::check_collinearity(model_lin) |>
    arrange(desc(VIF)) |>
    select(Term, VIF, VIF_CI_low, VIF_CI_high) |>
    performance::print_html()
Term VIF VIF_CI_low VIF_CI_high
lenke 29.92 24.23 37.00
regiao_cobb_inicial 22.13 17.95 27.34
dif_colete 6.47 5.33 7.92
correcao_colete 6.09 5.02 7.45
lordose_categ 3.64 3.05 4.41
cifose_categ 3.45 2.89 4.18
lordose_lombar 3.43 2.88 4.15
cifose_toracica 3.19 2.68 3.86
escoliometro_maior_10_graus 3.04 2.56 3.66
risser 2.42 2.06 2.90
idade 2.07 1.78 2.47
altura 1.64 1.44 1.94
sexo 1.40 1.25 1.65
flexibilidade 1.19 1.09 1.42
peso 1.16 1.07 1.39
Código
performance::check_collinearity(model_lin) |> plot()

2.1.2.2 Heteroscedasticidade

Verificação de heteroscedasticidade entre as variáveis independentes.

Código
performance::check_heteroscedasticity(model_lin)
Warning: Heteroscedasticity (non-constant error variance) detected (p = 0.017).
Código
performance::check_heteroscedasticity(model_lin) |> plot()

2.1.2.3 Outliers

Código
performance::check_outliers(model_lin)
OK: No outliers detected.
- Based on the following method and threshold: cook (0.921).
- For variable: (Whole model)
Código
performance::check_outliers(model_lin) |> plot()

2.1.2.4 Valores ajustados (diagnóstico)

Código
performance::check_predictions(model_lin) |> plot()

2.1.2.5 Resíduos simulados

Código
performance::check_residuals(model_lin)
OK: Simulated residuals appear as uniformly distributed (p = 0.753).
Código
performance::check_residuals(model_lin) |> plot()

2.1.3 Qualidade do ajuste do modelo

2.1.3.1 Coeficiente de determinação

\(R^2\).

Código
# Qualidade do ajuste
performance::r2(model_lin)
# R2 for Linear Regression
       R2: 0.369
  adj. R2: 0.293

2.1.3.2 Tamanho amostral efetivo (casos completos)

Código
tibble(n_modelo = nobs(model_lin))




2.2 Modelo de regressão logística

Este modelo estima associações ajustadas com a chance de melhora (MCID), definida como \(\Delta \le -5°\).

Código
model_bin <- glm(
    delta_cat ~
        # numericas
        idade +
        altura +
        peso +
        cifose_toracica +
        lordose_lombar +
        dif_colete +
        correcao_colete +

        # categoricas
        sexo +
        lenke +
        risser +
        flexibilidade +
        regiao_cobb_inicial +
        escoliometro_maior_10_graus +
        lenke_curvatura +
        risser_agrupado +
        cifose_categ +
        lordose_categ,
    data = df, family = "binomial"
)

2.2.1 Parâmetros do modelo

Parâmetros do modelo de regressão logística. Resposta em Odds Ratio.

Código
model_bin |>
    gtsummary::tbl_regression(exponentiate = TRUE, conf.level = 0.95) |>
    # gtsummary::add_global_p() |>
    gtsummary::bold_p() |>
    gtsummary::bold_labels() |>
    gtsummary::italicize_levels()
Characteristic OR 95% CI p-value
idade 0.94 0.71, 1.24 0.7
altura 0.90 0.01, 133 >0.9
peso 1.01 1.00, 1.04 0.3
cifose_toracica 1.04 0.99, 1.09 0.2
lordose_lombar 0.94 0.89, 0.99 0.012
dif_colete 0.88 0.78, 1.00 0.049
correcao_colete 1.02 0.98, 1.06 0.4
sexo


    feminino
    masculino 2.31 0.69, 7.94 0.2
lenke


    1
    2 0.12 0.02, 0.75 0.029
    3 0.07 0.01, 0.33 0.002
    4 0.04 0.00, 0.23 <0.001
    5 0.26 0.01, 5.89 0.4
    6 0.12 0.00, 2.24 0.2
risser


    0
    1 3.21 1.13, 9.48 0.031
    2 1.41 0.52, 3.85 0.5
    3 0.68 0.23, 1.94 0.5
    4 1.66 0.47, 6.02 0.4
flexibilidade


    flexivel
    rigido 0.45 0.22, 0.89 0.024
regiao_cobb_inicial


    lombar
    toracica 0.66 0.02, 8.61 0.8
    toracico_proximal 12.5 0.17, 1,130 0.2
escoliometro_maior_10_graus


    lombar
    nenhum 1.11 0.44, 2.74 0.8
    toracica 1.61 0.56, 4.70 0.4
    toracica_lombar 1.19 0.29, 4.96 0.8
lenke_curvatura


    curva em c
    curva em s


risser_agrupado


    final do crescimento
    inicio do crescimento


cifose_categ


    hipercifose
    hipocifose 1.65 0.14, 20.2 0.7
    normalidade 2.27 0.43, 12.2 0.3
lordose_categ


    hiperlordose
    hipolordose 0.16 0.04, 0.61 0.009
    normalidade 0.43 0.13, 1.39 0.2
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

2.2.2 Verificações de adequação do modelo

2.2.2.1 Colinearidade

Verificação de colinearidade entre as variáveis independentes (Variance Inflation Factor). Valores maiores que 10 indicam colinearidade entre as variáveis.

Código
performance::check_collinearity(model_bin) |>
    arrange(desc(VIF)) |>
    select(Term, VIF, VIF_CI_low, VIF_CI_high) |>
    performance::print_html()
Term VIF VIF_CI_low VIF_CI_high
lenke 33.43 27.06 41.35
regiao_cobb_inicial 25.05 20.31 30.96
dif_colete 4.76 3.95 5.80
correcao_colete 4.37 3.63 5.31
lordose_categ 3.78 3.16 4.59
lordose_lombar 3.44 2.88 4.16
cifose_categ 3.26 2.74 3.94
escoliometro_maior_10_graus 3.13 2.63 3.78
cifose_toracica 3.11 2.62 3.76
risser 2.52 2.14 3.02
idade 2.08 1.78 2.47
altura 1.65 1.44 1.95
sexo 1.51 1.33 1.78
flexibilidade 1.18 1.08 1.41
peso 1.17 1.07 1.40
Código
performance::check_collinearity(model_bin) |> plot()

2.2.2.2 Resíduos binarizados

Código
performance::binned_residuals(model_bin) |> plot()

2.2.2.3 Outliers

Código
performance::check_outliers(model_bin)
OK: No outliers detected.
- Based on the following method and threshold: cook (0.909).
- For variable: (Whole model)
Código
performance::check_outliers(model_bin) |> plot()

2.2.2.4 Valores ajustados (diagnóstico)

Código
performance::check_predictions(model_bin) |> plot()

2.2.2.5 Resíduos simulados

Código
performance::check_residuals(model_bin)
OK: Simulated residuals appear as uniformly distributed (p = 0.220).
Código
performance::check_residuals(model_bin) |> plot()

2.2.3 Qualidade do ajuste do modelo

Teste de Hosmer-Lemeshow para avaliar a qualidade do ajuste de modelos binomiais. P-valor < 0.05 indica que o modelo não ajusta bem os dados.

Código
performance::performance_hosmer(model_bin)
# Hosmer-Lemeshow Goodness-of-Fit Test

  Chi-squared: 13.960
           df:  8    
      p-value:  0.083

2.2.3.1 Coeficiente de determinação (Tjur)

Coeficiente de determinação de Tjur (\(R^2_{Tjur}\)). Teste específico para modelos binomiais.

Código
performance::r2(model_bin)
# R2 for Logistic Regression
  Tjur's R2: 0.334

2.2.3.2 Tamanho amostral efetivo (casos completos)

Código
mf_bin <- model.frame(model_bin)
bin_n <- nrow(mf_bin)
bin_events <- sum(mf_bin[[1]] == 1, na.rm = TRUE)
bin_nonevents <- sum(mf_bin[[1]] == 0, na.rm = TRUE)
tibble(
    n_modelo = bin_n,
    eventos_mcid = bin_events,
    nao_eventos = bin_nonevents
)




2.3 CART (árvore de decisão) para MCID (delta_cat)

2.3.1 Racional teórico

Um modelo CART (Classification and Regression Tree) é útil aqui porque:

  • Não linearidade: identifica automaticamente limiares (pontos de corte) em preditores contínuos.
  • Interações: a sequência de divisões (splits) representa interações condicionais sem precisar pre-especificar termos de interação.
  • Interpretabilidade: gera regras do tipo “se… então…”, úteis para comunicação clínica e geração de hipóteses.

Limitações importantes: árvores isoladas tendem a ter alta variância e podem superajustar. Por isso, esta seção usa validação cruzada (CV) e tuning/poda para controlar complexidade e reporta desempenho interno (sem validação externa).

2.3.2 Setup e preparação dos dados

Código
library(tidymodels)
library(rpart.plot)
library(vip)

set.seed(123)

# Desfecho como fator: "melhora" (evento) vs "nao_melhora"
df_cart <- df |>
    mutate(
        delta_cat_cart = factor(
            delta_cat,
            levels = c(0, 1),
            labels = c("nao_melhora", "melhora")
        )
    ) |>
    select(
        delta_cat_cart,
        # numericas
        idade,
        altura,
        peso,
        cifose_toracica,
        lordose_lombar,
        dif_colete,
        correcao_colete,

        # categoricas
        sexo,
        lenke,
        risser,
        flexibilidade,
        regiao_cobb_inicial,
        escoliometro_maior_10_graus,
        lenke_curvatura,
        risser_agrupado,
        cifose_categ,
        lordose_categ,
    )

# Verificar niveis (evento = "melhora", segundo nivel)
levels(df_cart$delta_cat_cart)
[1] "nao_melhora" "melhora"    

2.3.3 Recipe (pre-processamento)

O pre-processamento é feito dentro do resampling (evita leakage):

  • step_zv(): remove variáveis com variância zero
  • step_unknown() + step_novel(): robustez a níveis novos/raros durante CV
  • step_impute_median() + step_impute_mode(): imputação de missings
Código
cart_rec <- recipe(delta_cat_cart ~ ., data = df_cart) |>
    step_zv(all_predictors()) |>
    step_unknown(all_nominal_predictors()) |>
    step_novel(all_nominal_predictors()) |>
    step_impute_median(all_numeric_predictors()) |>
    step_impute_mode(all_nominal_predictors())

cart_rec

2.3.4 Modelo e hiperparâmetros

Hiperparâmetros a serem tunados:

  • cost_complexity (cp): controla poda/complexidade; valores maiores = árvores menores.
  • tree_depth: limita profundidade (proxy do grau de interações).
  • min_n: mínimo de observações em no terminal; evita regras baseadas em poucos pacientes.
Código
cart_spec <- decision_tree(
    cost_complexity = tune(),
    tree_depth = tune(),
    min_n = tune()
) |>
    set_engine("rpart") |>
    set_mode("classification")

cart_wf <- workflow() |>
    add_recipe(cart_rec) |>
    add_model(cart_spec)

cart_wf
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: decision_tree()

── Preprocessor ────────────────────────────────────────────────────────────────
5 Recipe Steps

• step_zv()
• step_unknown()
• step_novel()
• step_impute_median()
• step_impute_mode()

── Model ───────────────────────────────────────────────────────────────────────
Decision Tree Model Specification (classification)

Main Arguments:
  cost_complexity = tune()
  tree_depth = tune()
  min_n = tune()

Computational engine: rpart 

2.3.5 Validação cruzada e tuning

Usamos CV estratificada repetida (10-fold, 5 repetições) para estimar desempenho e selecionar hiperparâmetros.

Código
set.seed(123)
folds <- vfold_cv(df_cart, v = 10, repeats = 5, strata = delta_cat_cart)

# Metricas de interesse
metricas <- metric_set(roc_auc, accuracy, sens, spec)

# Grid de hiperparametros
set.seed(123)
cart_grid <- grid_latin_hypercube(
    cost_complexity(range = c(-4, -1)),
    tree_depth(range = c(3L, 6L)),
    min_n(range = c(5L, 30L)),
    size = 30
)

# Tuning
set.seed(123)
cart_tuned <- tune_grid(
    cart_wf,
    resamples = folds,
    grid = cart_grid,
    metrics = metricas,
    control = control_grid(save_pred = TRUE)
)

2.3.6 Melhores hiperparâmetros

Selecionamos pelo maior AUC ROC medio na CV:

Código
# Top 10 combinacoes por AUC
cart_tuned |>
    collect_metrics() |>
    filter(.metric == "roc_auc") |>
    arrange(desc(mean)) |>
    slice_head(n = 10)
Código
# Melhor combinacao
best_params <- cart_tuned |>
    select_best(metric = "roc_auc")

best_params

2.3.7 Desempenho final (CV)

Avaliamos o modelo finalizado com os melhores hiperparâmetros:

Código
cart_final_wf <- finalize_workflow(cart_wf, best_params)

set.seed(123)
cart_res <- fit_resamples(
    cart_final_wf,
    resamples = folds,
    metrics = metricas,
    control = control_resamples(save_pred = TRUE)
)

# Metricas agregadas (media e SD)
cart_metrics <- cart_res |>
    collect_metrics() |>
    select(.metric, mean, std_err, n) |>
    mutate(
        ci_lower = mean - 1.96 * std_err,
        ci_upper = mean + 1.96 * std_err
    )

cart_metrics

2.3.8 Métricas por fold (dispersão)

Visualização da incerteza nas estimativas de desempenho:

Código
cart_res |>
    collect_metrics(summarize = FALSE) |>
    ggplot(aes(x = .metric, y = .estimate, fill = .metric)) +
    geom_boxplot(alpha = 0.7, show.legend = FALSE) +
    geom_jitter(width = 0.1, alpha = 0.3, size = 1) +
    scale_fill_brewer(palette = "Set2") +
    labs(
        title = "Distribution of metrics by fold (CV repeated)",
        subtitle = "10-fold CV x 5 repetitions = 50 estimates per metric",
        x = "Metric",
        y = "Value"
    ) +
    theme_minimal(base_size = 12) +
    theme(
        plot.title = element_text(face = "bold"),
        axis.text.x = element_text(angle = 0)
    )

2.3.9 Matriz de confusão agregada (normalizada)

Código
# Predicoes agregadas de todos os folds
preds_all <- cart_res |>
    collect_predictions()

# Matriz de confusao
conf_mat_data <- preds_all |>
    conf_mat(truth = delta_cat_cart, estimate = .pred_class)

# Normalizar por linha (proporcao por classe verdadeira)
conf_mat_tbl <- conf_mat_data$table |>
    as.data.frame() |>
    group_by(Truth) |>
    mutate(
        prop = Freq / sum(Freq),
        label = scales::percent(prop, accuracy = 0.1)
    ) |>
    ungroup()

# Visualizacao heatmap normalizado
ggplot(conf_mat_tbl, aes(x = Prediction, y = Truth, fill = prop)) +
    geom_tile(color = "white", linewidth = 1) +
    geom_text(aes(label = label), size = 5, fontface = "bold") +
    scale_fill_gradient(
        low = "white",
        high = "steelblue",
        labels = scales::percent,
        name = "Proporcao"
    ) +
    scale_y_discrete(limits = rev) +
    labs(
        title = "Confucion Matrix",
        subtitle = "Proportion by true class (rows sum to 100%)",
        x = "Prediction",
        y = "True"
    ) +
    theme_minimal(base_size = 12) +
    theme(
        plot.title = element_text(face = "bold"),
        panel.grid = element_blank()
    )

2.3.10 Árvore final

Ajustamos o modelo final no dataset completo para visualização e interpretação:

Código
cart_fit_full <- fit(cart_final_wf, data = df_cart)
cart_rpart <- extract_fit_engine(cart_fit_full)

# Plot da arvore com rpart.plot
rpart.plot(
    cart_rpart,
    type = 4,
    extra = 104,
    under = TRUE,
    fallen.leaves = TRUE,
    roundint = FALSE,
    box.palette = "BuGn",
    shadow.col = "gray",
    main = "Decision Tree for MCID (delta <= -5 degrees)"
)

2.3.11 Importância de variáveis

Código
# Grafico de importancia com vip
cart_fit_full |>
    extract_fit_parsnip() |>
    vip(
        num_features = 15,
        geom = "col",
        aesthetics = list(fill = "steelblue", alpha = 0.8)
    ) +
    labs(
        title = "Variable Importance (CART)",
        subtitle = "Based on impurity reduction (Gini)",
        x = "Importance",
        y = "Variable"
    ) +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(face = "bold"))

2.3.12 Tabela de folhas (grupos de predicao)

Esta tabela mostra os nós terminais (folhas) da árvore, com as regras que definem cada grupo, o tamanho do grupo e a proporção de melhora:

Código
# Extrair informacoes dos nos terminais
frame <- cart_rpart$frame
leaves_idx <- which(frame$var == "<leaf>")

# Funcao para extrair regras de cada no
get_path_rules <- function(tree, node_id) {
    path <- rpart::path.rpart(tree, nodes = node_id, print.it = FALSE)
    rules <- path[[1]][-1] # remover "root"
    if (length(rules) == 0) {
        return("(raiz)")
    }
    paste(rules, collapse = " E ")
}

# Construir tabela de folhas
leaves_tbl <- tibble(
    no = as.integer(rownames(frame)[leaves_idx]),
    n = frame$n[leaves_idx],
    n_melhora = frame$n[leaves_idx] * frame$yval2[leaves_idx, 5],
    n_nao_melhora = frame$n[leaves_idx] * frame$yval2[leaves_idx, 4],
    prop_melhora = frame$yval2[leaves_idx, 5],
    predicao = ifelse(frame$yval[leaves_idx] == 2, "melhora", "nao_melhora")
) |>
    rowwise() |>
    mutate(regras = get_path_rules(cart_rpart, no)) |>
    ungroup() |>
    mutate(
        prop_melhora_pct = scales::percent(prop_melhora, accuracy = 0.1),
        grupo = paste0("Grupo ", row_number())
    ) |>
    select(grupo, predicao, n, prop_melhora_pct, regras) |>
    arrange(desc(predicao), desc(n))

leaves_tbl |>
    rename(
        Grupo = grupo,
        Predicao = predicao,
        N = n,
        `% Melhora` = prop_melhora_pct,
        Regras = regras
    ) |>
    knitr::kable(align = c("l", "l", "r", "r", "l"))
Grupo Predicao N % Melhora Regras
Grupo 1 nao_melhora 147 27.2% dif_colete>=-19.5 E lenke=2,3,4,6
Grupo 2 nao_melhora 8 12.5% dif_colete>=-19.5 E lenke=1,5 E dif_colete>=-15.5 E flexibilidade=rigido
Grupo 5 melhora 75 77.3% dif_colete< -19.5
Grupo 3 melhora 12 66.7% dif_colete>=-19.5 E lenke=1,5 E dif_colete>=-15.5 E flexibilidade=flexivel
Grupo 4 melhora 10 100.0% dif_colete>=-19.5 E lenke=1,5 E dif_colete< -15.5

2.3.13 Curva ROC

Código
# Curva ROC a partir das predicoes CV
roc_data <- preds_all |>
    roc_curve(truth = delta_cat_cart, .pred_melhora, event_level = "second")

# AUC para o titulo
auc_val <- preds_all |>
    roc_auc(truth = delta_cat_cart, .pred_melhora, event_level = "second") |>
    pull(.estimate)

ggplot(roc_data, aes(x = 1 - specificity, y = sensitivity)) +
    geom_path(linewidth = 1.2, color = "steelblue") +
    geom_abline(linetype = "dashed", color = "gray50") +
    geom_ribbon(aes(ymin = 0, ymax = sensitivity), alpha = 0.1, fill = "steelblue") +
    annotate(
        "text",
        x = 0.6, y = 0.3,
        label = paste0("AUC = ", round(auc_val, 3)),
        size = 5, fontface = "bold"
    ) +
    labs(
        title = "ROC Curve (CV)",
        subtitle = "Based on aggregated predictions from all folds",
        x = "1 - Specificity (False Positive Rate)",
        y = "Sensitivity (True Positive Rate)"
    ) +
    coord_equal() +
    theme_minimal(base_size = 12) +
    theme(plot.title = element_text(face = "bold"))

2.3.14 Interpretação

Como ler a árvore:

  • Cada nó interno representa uma condição (split) sobre uma variável.
  • Cada caminho da raiz ate um nó terminal representa uma regra de classificação.
  • A sequencia de splits indica interações condicionais: um split so importa dentro da região definida por splits anteriores.
  • Os limiares (pontos de corte) sugerem não-linearidades nos efeitos.

Como ler a importância:

  • Variáveis com maior importância contribuem mais para a redução de impureza (Gini) nos splits.
  • Isso não implica causalidade, apenas relevância preditiva no contexto do modelo.

Limitações:

  • Desempenho reportado e interno (CV), sem validação externa.
  • Árvores são instáveis: pequenas mudanças nos dados podem alterar a estrutura.
  • Usar como ferramenta exploratória complementar aos modelos inferenciais (linear/logístico).